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Abstract 

In this work we present a detailed study of the Fermion Monte Carlo algorithm (FMC), a recently 
proposed stochastic method for calculating fermionic ground-state energies. A proof that the FMC 
method is an exact method is given. In this work the stability of the method is related to the 
difference between the lowest (bosonic-type) eigenvalue of the FMC diffusion operator and the 
exact fermi energy. It is shown that within a FMC framework the lowest eigenvalue of the new 
diffusion operator is no longer the bosonic ground-state eigenvalue as in standard exact Diffusion 
Monte Carlo (DMC) schemes but a modified value which is strictly greater. Accordingly, FMC 
can be viewed as an exact DMC method built from a correlated diffusion process having a reduced 
Bose-Fermi gap. As a consequence, the FMC method is more stable than any transient method (or 
nodal release-type approaches). It is shown that the most recent ingredient of the FMC approach 
[M.H. Kalos and F. Pederiva, Phys. Rev. Lett. 85, 3547 (2000)], namely the introduction of 
non-symmetric guiding functions, does not necessarily improve the stability of the algorithm. We 
argue that the stability observed with such guiding functions is in general a finite-size population 
effect disappearing for a very large population of walkers. The counterpart of this stability is 
a control population error which is different in nature from the standard Diffusion Monte Carlo 
algorithm and which is at the origin of an uncontrolled approximation in FMC. We illustrate the 
various ideas presented in this work with calculations performed on a very simple model having 
only nine states but a full "sign problem" . Already for this toy model it is clearly seen that FMC 
calculations are inherently uncontrolled. 

PACS numbers: 02.70Ss, 05.30.Fk 
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I. INTRODUCTION 



In theory quantum Monte Carlo (QMC) techniques[l| are capable of giving an exact 
estimate of the energy with an evaluation of the error: the statistical error. Unfortunately, 
such an ideal situation is not realized in practice. Exact results with a controlled finite 
statistical error are only achieved for bosonic systems. For fermionic systems, we do not 
have at our disposal an algorithm which is both exact and stable (statistical fluctuations 
going to zero in the large simulation time regime). This well-known problem is usually 
referred to as the "sign problem". The usual solution to cope with this difficulty consists 
in defining a stable algorithm based on an uncontrolled approximation, the so-called Fixed- 
Node approximation. {2], s], li, 1^, 0, In practice, the fixed-node error on the energy is 
small when one uses good trial wavefunctions and, thus, QMC methods can be considered 
today as reference methods to compute groundstate energies as shown by a large variety of 
applications fl, S S Q, [n , Q, [13I . However, the accuracy of the results is never known from 
the calculation, it is known only a posteriori, for example by comparison with experimental 
data. Exact methods, which are including the 

nodal release method[lJ], have been applied with success only to very specific models (small 
or homogeneous systems) for which the sign instability is not too severe (small Bose-Fermi 
energy gap). 

Recently, Kalos and coworkers 18| have proposed a method presented as curing the sign 
problem, the so-called Fermion Monte Carlo method (FMC). This work makes use of two 
previously introduced ingredients, a cancellation process between "positive" and "negative" 
walkers introduced by Arnow et al. [l^ and a modified process correlating explicitly the 



dynamics of the walkers of different signs. 



201 ] The new feature introduced in 18| is the 



introduction of non-symmetric guiding functions. The method has been tested on various 
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2l| 



simple systems including free fermions and interacting systems such as the ^He fiuid 
The results are found to be compatible with the assumption that the method is stable and 
not biased. However, this conclusion is not clear at all because of the presence of large 
error bars. The purpose of this work is to present a detailed analysis of the algorithm and 
a definitive answer to this assumption. 

The content of this paper is as follows. In section II we give a brief presentation of the 
"sign problem" . The sign instability in an exact DMC approach comes from the blowing up 
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in imaginary time of the undesirable bosonic component associated with the lowest math- 
ematical eigenstate of the Hamiltonian (a wavefunction which is positive and symmetric 
with respect to the exchange of particles). The fluctuations of the transient energy esti- 
mator grows like e^^^o-^o)^ where Eq is the fermionic ground-state energy (here and in 
what follows, the superscript F stands for "Fermionic") and Eq , the lowest mathematical 
eigenvalue of the Hamiltonian (the superscript "B" standing for Bosonic). In Sec. Ill we 
briefly recall the main elements of the standard DMC method, the basic stochastic algorithm 
simulating the imaginary-time evolution of the Hamiltonian. In section IV, we describe the 
FMC method as a generalization of the DMC method and we show that FMC is an exact 
method, that is, no systematic bias is introduced. This section is made of two parts. In 
a first part we introduce the notion of positive and negative walkers to represent a signed 
wavefunction in DMC. This will help us to view the FMC method as a generalization of the 
DMC method, the FMC introducing two important modifications with respect to DMC: a 
correlated dynamics for positive and negative walkers and a cancellation process for such 
pairs whenever they meet. 

In section V we study the stability of the algorithm. We prove that the FMC method 
is in general not stable, the fiuctuations of the transient estimator of the energy growing 
exponentially like e^'^^^-^o) where E^ is the lowest bosonic-like eigenvalue of the generalized 
diffusion process operator associated with FMC whose expression is given explicitly. It 
is shown that FMC is more stable than the standard nodal release DMC method because 
Eq > Eq . In section VI we illustrate our theoretical results using a toy model, a "minimal" 
quantum system having a genuine sign problem (two coupled oscillators on a finite lattice). 
The different aspects of the FMC method are highlighted in this application. In the case 
of non-symmetric wavefunctions introduced recently in Ref. [isl it is shown that, in contrast 
with DMC where the population control error decays linearly as a function of the population 
size, the FMC decay displays a much more slower power law. As a consequence, one needs a 
very large population of walkers to remove the control population error and to observe the 
instability of FMC. Let us emphasize that observing such a subtle finite population effect 
for a genuine many-fermion system is actually impossible. Here, to illustrate this important 
point numerically, we have been led to consider a very simple system having only a few states. 
For this system, a large population of walkers -eventually much larger than the dimension 
of the quantum Hilbert space itself- can be considered. The results obtained in that regime 



4 



confirm our theoretical findings, in particular the fact that the stability of the algorithm 



presented elsewhere 



18| is only apparent. As an important conclusion, we emphasize that 



the FMC control population error is an uncontrolled approximation for realistic fermion 
systems. 



II. FERMION INSTABILITY IN QUANTUM MONTE CARLO 

We consider a Schrodinger operator for a system of fermions, 

H = -iv^ + ViR) (1) 

where we note R = {ri . . .r^) the 3A^ coordinates of the N particles in the three dimensional 
space. In this expression, the first term is the kinetic energy (V^ is the Laplacian operator in 
the space of the 3A^ coordinates). The second term, V, is the potential. DMC techniques are 
based upon the evolution of the Hamiltonian in imaginary time. We express this evolution 
using the spectral decomposition: 

e-*(^-^-) = Ee^*^'''"''"M0.)(0.l (2) 

i 

where 0j are the (normalized) eigenfunctions of H, Ei are the corresponding eigenvalues and 
Et is a so-called reference energy. The fundamental property of operator ([2]) is to filter out 
the lowest eigenstate 0o- To understand how it works, we consider the time evolution of a 
wavefunction /o(R) 

I /,) ^ e-*(^-^-) I /o) (3) 
and calculate the overlap with an eigenstate 0j 

(0. |/*) = e-*(^'-^-)(0. l/o). (4) 

From this expression it is easily seen that the component on the lowest eigenstate 0o is 
growing exponentially faster than the higher components. In DMC methods the function 
is generated using random walks. For t large enough the lowest eigenstate 

ft ~ e-*(^°-^-)0o (5) 

is produced. This asymptotic behaviour makes DMC an accurate method for computing 
the properties of in particular the energy Eq. Unfortunately, for fermionic systems the 
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physical groundstate 0^, which is antisymmetric with respect to the exchange of particles, 
is not (pQ, but some "mathematically" excited state because 0o is positive and symmetric for 
a Schrodinger operator (bosonic ground-state). For the sake of clarity, we shall denote from 
now on this bosonic ground-state as (f)Q and its energy, Eq. The necessity of extracting an 
exponentially small component to evaluate Eq or any fermionic property is at the origin of 
the fermion instability in quantum Monte Carlo. 

We now give a quantitative analysis of this instability. Since the asymptotic behaviour 
of ft is not useful to compute the fermionic e nerg y, we consider the transient behaviour of 
the evolution of ft- Basically, exact methods [14] filter out (pQ by projecting the transient 
behaviour of ft on the antisymmetric space. Introducing the antisymmetrization operator 
A, one obtains the physical groundstate for t large enough 

A I /,) ~ e-*(<-^-)< (6) 

since the components of ft over the higher antisjnnmetric eigenstates are decreasing expo- 
nentially with respect to the component on 0^. In practice, the fermionic energy Eq is 
calculated using an antisymmetric function ipT and using the fact that, at large enough t, 
one has 





H 


\ft) 




T 


\ft) 



(7) 



In the long-time regime the stochastic estimation of the R.H.S. of Eq.Q is unstable. In 
essence, the signal, the antisymmetric component of ft-, mettons des virgules.... decreases 
exponentially fast with respect to ft-, Eq.([5]). The signal-over- noise ratio behaves like 
^-t(E^-E^) QT^f^^ thus, an exponential growth of the fluctuations of the DMC estimator, 
Eq-dZD, appears. 

Now, let us give a more quantitative analysis by writing this estimator and evaluating the 
variance. In a standard diffusion Monte Carlo calculation, the time-dependent distribution 
ft is generated by a random walk over a population of walkers {Rj}. Formally ft is given 
in the calculation as an average at time t over Dirac functions centered on the walkers and 
weighted by some positive function ipc 

/.(R)-^(i:^(R-R<)> (8) 

In this expression, (...) denotes the average over populations of walkers {Rj} obtained at 
the given time t. The function ipc is usually called the importance or guiding function. 
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Replacing ft in ([7]) by its expression ([8]), the estimator of the energy reads for large enough 
t 

In practice, both numerator and denominator are computed as an average over the Ns 
walkers produced by the algorithm at time t. As a consequence, the energy is obtained as 

Ns ^«=1 V'G 

The ratio ^ is an estimator of the energy for A^^ large enough, when the numerator and 
denominator have small fluctuations aroung their average. Now, let us evaluate the fluctu- 
ations of this ratio in the limit of a large population Ns 

^ Kv) = W • ^ ^ 

Here, we have used the hypothesis that the fluctuations of the denominator and the numer- 
ator are very small. Using the fact that M and V are statistical averages over independent 
random variables with the same distribution one obtains 

1 



■<Pg ^ 



) 



We can replace these averages by integrals over the distribution ftipc, Eq-dH]) 



(12) 



v)-Ns (V^T^n^ at/ ^^^^ 



which confirms quantitatively that the statistical error grows exponentially in time. Let 
us emphasize that this problem is particularly severe because the Bose-Fermi energy gap, 
Ab~f = —Eq , usually grows faster than linearly as a function of the number of particles. 
In practice, one has to find a trade-off between the systematic error coming from short 
projection times t and the large fluctuations arising at large projection times. 



III. THE DIFFUSION MONTE CARLO METHOD 



The FMC method is a generalization of the well-known DMC method. Presenting this 
algorithm is a useful preparation for the next section about Fermion Monte Carlo. The 
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DMC method generates the function ft following the imaginary time dymamics of H 

f\ ^ e-*(^-^-)/o (14) 

where /o is a positive function. The imaginary time dynamics is produced by iterating many 
times the short-time Green function ^-'^iH-Er) ^^gj^g r is a small time step. The distribution 
ft> at the time t' = t + t is then obtained from ft as follows 

f, = e-{^-^-)/, (15) 

where the density ft is sampled by the population of walkers {Rj}. Using the Dirac ket 
notation, ft given by Eq.([H]) is rewritten as 

/.^^(eir.)). (16) 

where ipG is some positive function, the so-called guiding function. Let us show how the 
density ft' is generated from the distribution (fTB|) . Replacing in f[T^ the function ft by the 
R.H.S. of (USD one has 

ft' = e-^^^'-^-^^^ElR^)) (17) 
= (Ee-^"-^^^^|R.)) (18) 
^ ^E^^'IR^)^ (19) 



where we have introduced the operator 



L = -^^{H-Et)^. (20) 



For a Schrodinger Hamiltonian ([T]) the operator L takes the form 



L = -^g{H - El)^ - {El - Et) (21) 

Wo 

= lv^-V[b.] (22) 

-{El-Et) (23) 



where we have introduced the so-called drift vector 

b = — — 



(24) 



and the local energy of the guiding function ipc, 

^ (25) 

The operator L is the sum of the so-called Fokker Planck operator (!22|) and a local operator 
fl23|l . Using this decomposition, the vector e^^ \ Rj) appearing in the average f|T9l) can be 
rewritten for small enough time step r as follows 

I R.) = e-[^^'-^[b-]] (26) 

^^-r{E^-ET) I (27) 

The action of e^^ on | Rj) is sampled as follows. The short-time dynamics of the Fokker 
Planck operator fl26|) is performed by the way of a Langevin process, 

Rr = R^^ + b>tT + ^r^^ (28) 

where /i runs over the 3A^ coordinates (three space coordinate for each fermion), and T^f are 
independent gaussian random variables centered and normalized 

(r^fr/n = (29) 

The averaged Langevin process (!28|) is equivalent to apply the short-time dynamics of the 
Fokker Planck operator ( l26i) : 

R^)\ = e^[iv^-^N (30) 



The factor 



= e--(^^(i^i)-^^) (31) 



being a normalization term, called the branching term. The new walker R^ is duplicated 
(branched) a number of times equal to Wi in average. This process is a birth-death pro- 
cess since some walkers can be duplicated and some can be removed. The pop ulation of 
walkers fluctuating, one has to resort to control population techniques. [22|, |23|, l2J] With 



these two processes, diffusion and branching a new population of walkers {R^} is produced, 
representing in average the desired result : 

(EiR:)\ = (Ee^'iR^)\ (32) 



From (fT9|) and (132|) one can see that the distribution ff is sampled by the new population 
of walkers {Ri} according to f|T6|) 



(33) 



In summary, by iterating these two simple operations, namely the Langevin and branch- 
ing processes, the DMC method allows to simulate the imaginary time dynamics of the 
Hamiltonian, thus producing a sample of ft, Eq. flTBl) . Various properties of the system can 
be computed from this sample, e.g. ground-state bosonic energies, M, y, 
energies. |16l. I17| and various observables. [26|, 



27 



28 



291, lad. 
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251 ] excited-state 



For the vast majority of the DMC simulations on fermionic systems, only an approxima- 
tion of the exact fermionic ground-state energy is computed, namely the so-called Fixed- 



Node energy. 



fly. 
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35( 1 In a fixed-node DMC calculation the guiding function is 



chosen as ipc = |^r| where ipT is some fermionic antisymmetric trial wavefunction. With 
this choice, the guiding function vanishes at the nodes (zeroes) of the trial wavefunction and 
the drift vector diverges. As a consequence, the walkers cannot cross the nodes of ipT and 
are confined within the nodal regions of the configuration space. It can be shown that the 
resulting DMC stationary state is the best variational solution having the same nodes as 

rom the R.H.S of 



or ( jTOjl is an 
34| Note that, in practice, Eq^ is 



ipT- In other words, the "Fixed-Node" energy obtained 
upper bound of the exact fermi energy, Eq'^ > -Ejf . jssi 
in general a good approximation of the true energy, [s], IssI 

In the present work, we are considering "exact" DMC approaches for which the exact 
fermionic energy calculated from expression iQ is searched for. As discussed in the previous 
section, such exact DMC calculations are fundamentally unstable. A famous example of an 
exact DMC approach is the nodal release method of Ceperley and Alder, jl^ . Basically, 
nodal release methods are standard DMC methods where the fixed-node distribution (sam- 
pled with a standard Fixed-Node DMC) is chosen as initial distribution /q. In exact methods 
the guiding function ipc is strictly positive everywhere, so that the walkers can cross the 
nodes of the trial wavefunction. Exact fermi methods are efficient in practice only when the 
convergence of the estimator is fast enough, that is, when it occurs before the blowing up of 
fluctuations, Eq. (fT3l) . In practice, two conditions are to be satisfied. First, the fixed-node 
wavefunction must be already close enough to the exact solution 0^. For this reason 
the choice of the trial function ipx (quality of the nodes of ipx) is crucial. Second, the Bose- 
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Fermi gap, Ab-f = Eq — E^, which drives the asymptotic behaviour of the fluctuations, 
Eq. f|T3|) . must be small. The quantity — Eq depends only on the Hamiltonian at hand; 
there is no freedom in the nodal release method to modify the asymptotic behaviour. We 
will deflne in the next section the Fermion Monte Carlo method (FMC) as a generalization 
of the DMC method and show in sections V and VI that, in contrast with the nodal release 
method, the FMC method can improve substantially the asymptotic behaviour, Eq. (fT3|) . 

IV. THE FMC METHOD 

A. Preliminary: Introducing positive and negative walkers in DMC 

In Fermion Monte Carlo a dynamics on a signed function ft is performed. In what follows, 
we show that DMC can be easily generalized to the case of a signed distribution ft and, 
thus, FMC can be viewed as a simple generalization of DMC. If ft carries a sign, it can be 
written as the difference of two positive functions 

ft = f^ - ft- (34) 
both satisfying the following equations of evolution 

/+ = e-*(^-^-)/o+ (35) 
ft- = e-*(^-^-Vo" (36) 

To sample these expressions, two independent DMC calculations can be carried out. The 
positive part ft' is then sampled by a population of walkers {Ri^} (called "positive" walk- 
ers). The distribution ft~ is related to {Ri"*"} as in Eq. flT^ 

(37) 

and the negative part is similarly sampled by a population of "negative" walkers {Ri^} 

/.-^(EIHd). (38) 

Note that we consider here the general case where the guiding functions associated with the 
positive and negative walkers, ipQ and ipQ, are different. Finally, the dynamics of positive 
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and negative walkers is described by the DMC-like diffusion operators 



wliere tlie drift vectors are given by 



-^l;^{H - Et)-L (39) 

iv^-V[b±] (40) 
-{Et - Et) (41) 



and tlie local energies of the guiding functions ij^Q by 

Et.f. (43) 

In actual calculations, two Langevin dynamics on the positive and negative walkers are 
performed 

= i?/± + Ifr + ^/^r/f ^ (44) 
and positive and negative walkers are branched according to their respective weight 

^±(R±) = e-(^^*(i^?)-^-). (45) 
B. The detailed rules of FMC 

In a few words, the FMC method is similar to a DMC method on a signed function, 
except that the positive and negative walkers are correlated and can annihilate whenever 
they meet. The Langevin processes are correlated in such a way that positive walkers and 
negative walkers meet as much as possible. A cancellation procedure is then performed 
when a positive walker and a negative walker meet. We will see in the next section that 
this cancellation procedure is at the origin of an improved stability of the algorithm. In this 
section we give a complete description of the algorithm and show that this algorithm does 
not introduce any systematic bias. 

In FMC the guiding functions V'g V'g ^^t arbitrary, they are related under any 
permutation P of two particles as follows 

^+(R)=^g(^R)- (46) 
12 



Various choices are possible for the guiding functions. Here, we consider the form proposed 
by Kalos and Pederiva,[l8|] namely 

= V^Vj+^±c^r (47) 

where ips is a symmetric (Bose-like) function, ipx an antisymmetric trial wavefunction, and 
c some positive mixing parameter allowing to introduce some antisymmetric component into 

Having in mind this choice for the guiding functions we will show how the two DMC 
processes over the two populations of walkers {R/^} and {R^} can be replaced by a diffusion 
process over a population of pairs of walkers {(Rj^,Rj^)}. We first show how the DMC 
process can be modified to maintain as many positive and negative walkers during the 
simulation. In the DMC dynamics, the branching terms associated with the positive and 
negative walkers are in general different. As a consequence, the number of positive walkers 
A^^, can be different from the number of negative walkers A^^ . At time t the DMC density 
ft reads 

^-(S^^"^'>-S^"^^"V- '''' 

This formula is obtained by replacing in equation fl51|) the expressions f l57|) and fl55]) for 
and respectively. If and A^^ are different, we will replace ft (1481) by a new function 
Qt sampled with an equal number of positive and negative walkers. Such an operation does 
not introduce any bias if the antisymmetric components of the future evolution of ft and gt 
are identical. Indeed, only the antisymmetric component of ft contributes to the estimator 
of the energy, Eq. ([7]). At time t' > t the two densities are given by 

f, = e-(*'-*)(^-^-)/, (49) 
9, = e-(*'-*)(^-^-)<7*. (50) 

Let us write that the antisymmetric components of ft' and Qf must be equal 

Using the fact that the antisymmetrisation operator A commutes with the evolution operator 
and regrouping all the terms one finally finds 



e 



<t'-t){H- 



""^^Aift - gt) = 0. (52) 
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This condition is satisfied whenever A{ft — gt) = 0. The important conclusion is that one 
can replace ft by any function gt such that the difference gt — ft is orthogonal to the space 
of antisymmetric functions. Let us now show how this property can be used to impose a 
common number of walkers in the positive and negative populations. 

Let us consider the case where there are more positive walkers than negative walkers, 
Ng > Ng . In this case one can substract from ft, Eq. (H8|) . the following vector 

^|R.->+P^|R« (53) 

where is a positive walker and P is a two-particle permutation. Such an operation is 
allowed since the application of the antisymmetrizer to the vector fl53p gives zero [a direct 
consequence of .4(1 + P) = 0]. Now, because of Eq. fH^ the vector can also be written 

as 

4 I + ^/ 1 R?). (54) 

Substracting this vector from ft removes the contribution -V | R^) from ( HHi) and adds the 
contribution —-j^P \ R^^) to (HSj) . In other words, the positive walker R^ has been removed 
and the negative walker 

Rr = PR+ (55) 

has been created. Similarly, one can remove a negative walker R^" and create a positive 
walker 

R+ = PR-. (56) 

This possibility of transfering one walker from one population to the other one allows to 
keep an identical number of walkers for the two populations at each step. Now, thanks to 
this possibility, it is possible to interpret the two populations consisting of the Ns positive 
walkers R^ and the Ns negative walkers R" {i G [l..A^s]) as an unique population of Ns 
pairs of walkers {(Rj^,Rj")}. Following this interpretation, the density ft (l48l) can be then 
rewritten as an average over a population of pairs of walkers 

and the energy can be computed as a ratio of averages performed on the population of pairs 



E.(^(R.+) - ^(Rr)) 

= ^-^ ^, (58) 
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where Eq.(I7]) has been rewritten by replacing ft using Eq. (l57Il . Now, everything is in order 
to detail the short-time dynamics of FMC The FMC dynamics consists of three steps 
(Langevin, branching, and cancellation steps): 

(i) Langevin step The Langevin processes (jH]) are simulated as in DMC, except that the 
gaussian random variables of the positive and negative walkers are no longer independent. 
The positive walker and the negative walker are moved according to Eq.( 



R/^^ = + b'^^r + v^r/f ^ (59) 
where T^f ^ are gaussian centered random variables verifying 

ivtvt) = (60) 

Such a move insures that the density of positive and negative walkers obey the Fokker Planck 
equation: 



However, the gaussian random variables are no more independent, they are correlated 
within a pair 

cr ^ {v^vn ^ 0. (62) 

Differents ways of correlating positive and neg ative walkers can be considered. We shall 



employ here the approach used in Refs. 18|, l20| which consists in obtaining the vector t]. 



n 1 



representing the 3A^ coordinates rj^ , from fj^ by reflexion with respect to the hyperplane 
perpendicular to the vector R^^ — Rj". 

r^^- = rfi" - ^ %^^^fy (^t - RD- (63) 

This relation between the gaussian random variables makes the move deterministic along 
the direction R^*" — Rj" . Such a construction insures that the walkers within a pair will meet 
each other in a finite time (even in large-dimensional spaces) . This aspect will be illustrated 
numerically in the last section. Formally, the two correlated Langevin processes can be seen 
as one Langevin process in the space of pairs of walkers 

(R+',Rr') = (R+,Rr) + (b+(R+),b-(Rr))r + Vf(r^,rr) (64) 
where ff^ and rfl are related via (|63|1 . 
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(ii) Branching step As we have already noticed, the branching of the negative and the 
positive walker are different within a pair: 



w- 



Taking into account their respective weights, the two walkers of a pair give the following 
contribution to the density ft 

If for example w'l > w~ , this vector can be written as 



Ilt)-w7—\R7). (66) 



I R,+) - I Rr) (67) 
+ iwt-w^)^\^t) (68) 

rG 

This density is the sum of two contributions. The first contribution, Eq. fl67|) . comes from 
a pair of walkers (R^,R.7) and carries the weight . The second, Eq. (!68|) . comes from 
a single positive walker R^*" and carries the weight w^' — . This single walker R^ can 
be replaced by a pair as follows. First, this single walker can be replaced by two positive 
walkers R^ with half of the weight, ^iw^ — w~). One of these two positive walkers carrying 
half of the weight can be transfered to the population of negative walkers by exchanging two 
particles. Finally, this single walker R^*" can be replaced by a pair (Rf, PHt) carrying the 
weight ^{w^ — w^). The resulting process just described is a branching of the pair (R^, R^) 
with the weight w~ and the creation of the pair (R^,PRj^) with the weight ^{wl' — w^). 
Of course, if one has wf < w~, then, the pair (R^, R^") is branched with the weight wf 
and the pair (PRj~, Rj") is created with the weight \{w^ — w^). 

Both cases, < wf or > wf, can be summarized as follows. The pair of walkers 
(R^, R~) is branched with the weight 

mm{wt, w;) = e-(--(^J'^r)-^^) (69) 

and the pairs (R^,PR/') and (PR~,R^) are created with their respective weights 

^{wt - min(«;+, ^r)) = ^[Et - max(E+, EZ)] + 0{r') (70) 



and 



\{w- - min(«;+, wi)) = ^lEt - max(i?+, EZ)] + Oir') (71) 
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(in) Cancellation step The third step is a cancellation procedure performed whenever a 
positive and negative walker meet. When Ri^ = Ri~, the contribution of the pair to the 
density can be simplified as follows 



(72) 



If the term in brackets is positive, this contribution comes from one single positive walker 



Hf with multiplicity 



i-#(R^ 



One can transform this single walker into a pair of 



positive and negative walker (R^jPRj*") with the new multiplicity 

1 



(73) 



If the term in brackets is negative, the pair {PRl, Rl) is drawn and the multiplicity is given 
by 



R 



(74) 



This is a cancellation procedure because a pair (R^, R^) with a multiplicity 1 has been 
transformed into a pair with a multiplicity smaller than one. Note that, when iJjq = ipQ, the 
multiphcities ( 1731) or (17^ reduce both to zero. In other words, there is a total cancellation of 
the pair whenever the walkers meet. As we shall see in the next section, the cancellation step 
is at the origin of the improved stability. The basic reason is that this procedure removes 
pairs which do not contribute to the signal but only to the statistical noise. A rigorous 
analysis of this point is provided in the next section. 

V. STABILITY OF THE FMC METHOD 
A. Criterium for stability 



We have just seen that the FMC method is a generalization of the DMC approach and 
we have shown that FMC preserves the evolution of the antisymmetric component of the 
sampled density. Now, having shown that FMC is an exact method, it is necessary to study 
the stability of the method. For that purpose, we consider the estimator of the energy, Eq. 
(!58|) . in the large-time regime 



^0 = — - 



"Pa ^ 



(75) 
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^2 



where Ns is the population size at time t. In the same way as done for the estimator (fTOl) . 
we can evaluate the variance of Eq by supposing that Ns is large enough so that both 
numerator and denominator have small fluctuations around their average 

Let us begin with the denominator. Using identity ( 1571) . the average of the random variable 
V defined in (175|) is nothing but 

(1^) = ^(V^T I /*). (77) 

Replacing ft by its asymptotic behaviour, Eq.([6]), we finally find that the denominator of 
fl76|) behaves for large t as 

(P)^ = i,e-^*(<-^-)(^^|0^)2. (78) 

Now, let us compute the numerator of ( 1761) . This numerator can be written as the variance 
of a sum of random variables defined over the pairs of walkers 



{(Af-E^vr) 



E.r(R+„R: 



1 2\ 



Ns 

where we have introduced the function FfR"^, R^) 



(79) 



r(R.. R-) . - (£^|0*i(R-). (80) 

Using the fact that the pairs of walkers have the same distribution and supposing that they 
are independent we finally find 

((AT - E^Vf) = i-a^ ( r(R+, R-)) . (81) 

If the pairs of walkers are not independent this expression is only modified by a correlation 
factor independent on the time and the population size Ns provided that Ns and t are large 
enough. Finally, up to a multiplicative factor, the variance of the energy estimator has the 
following asymptotic behaviour 

a^{^) oc -^^C{t) (82) 
^V^ Ns{t) ^ ^ 

where the coefficient C{t) is given by 

C{t) = Ns{tfe-^'^''o-ET) (83) 
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and where Ns{t), the number of pairs, depends on time t due to the birth-death process. 
In expression fl83l) we note that the behaviour of fl76l) at large times t is related to the 
value of Eq and the asymptotic behaviour of the number of pairs of walkers Ns{t). We can 
already understand physically the interest of the cancellation process: this process limits 
the growth of Ns{t) the number of pairs of walkers, thus limiting the growth of C{t) (155]) 
and the variance (182|) . Let us now precise this criterium more rigorously by evaluating 
the asymptotic behaviour of Ns{t). For that purpose we introduce the density of pairs 
Utijlt, R-r); this density obeys a diffusion equation we write down 



-(Dfmc - ET)Ut (84) 



dt 

where we have introduced the diffusion operator — (-Dfmc — Et). We will give its expression 
later; for the present purpose we just need the asymptotic behaviour of lit given by 

Tit = e-*(^o"-^-)n5 (85) 

where lis is the stationary density of the process, namely the lowest eigenstate of the 
operator -Dfmc, and E^ is the corresponding eigenvalue. The number of pairs Nsif) behaves 
as the normalization of lit, and consequently grows like e"**^^?"^^^ Note that in practice one 
adjusts the reference energy Et to Eq during the simulation to keep a constant population of 
average size Ns along the dynamics. Such a procedure is referred to as a control population 



technique 



241 ] and will be discussed later. Finally, the asymptotic behaviour of the variance 



of the FMC estimator of the energy is 

^V' Ns 

This expression is analogous to fll3p except that the lowest energy of the Hamiltonian op- 
erator H has been replaced by the lowest energy of the operator -Dfmc- In conclusion the 
stability of the algorithm is related to the lowest eigenvalue of the FMC diffusion operator, 
Eq . It is clear from (186|) that the higher this eigenvalue is, the more stable the simulation 
will be. In the next section, we will discuss the allowed values of Eq . This will prove that 
FMC is not a stable method in general, but is more stable than any standard transient 
method. 
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B. Stability of the Fermion Monte Carlo algorithm 



In this section we prove that the lowest eigenvalue of the FMC operator, Eq , has the 
following upper and lower bounds 



^From the expression of the variance, Eq. fl86|) . one can easily understand the meaning of 
these two inequalities. The first inequality indicates that FMC is not a stable method, the 
stability being achieved only in the limit Eq = Eq . Note that, even for very simple systems, 
this stability is in general not obtained. This important point will be illustrated in the next 
section. The second inequality shows that FMC is more stable than any standard transient 
DMC method (nodal release method). Indeed, the exponent associated with the explosion 
of fluctuations, Eq. (!86|) . is smaller than in the standard case, Eq.([T3l). Before giving a 
mathematical proof of these two inequalities, let us first present some intuitive arguments 
in their favor. The first inequality, Eq. (1871) . takes its origin in the fact that the signal -the 
antisymmetric component of ft, Eq.([H]) is extracted from the population of pairs of walkers, 
Eq. llFri) . and, consequently, cannot grow faster than the population of pairs itself, Eq. (155i) . 
The second inequality can be understood as follows. Without the cancellation process, the 
FMC method reduces to two correlated DMC algorithms. The number of walkers grows as 
in a standard DMC, namely ~ e^i^o-^T)t _ ^pj^g cancellation process obviously reduces the 
growth of the population of walkers, Eq. fISS]) . and, thus, we should expect that Eq > Eq. 
Now, let us give some more rigorous proofs. For that purpose we compare the Fermion 
Monte Carlo operators with and without cancellation process. Without the cancellation 
process the FMC diffusion operator reads 



E^ < < 



(87) 
(88) 








D-Et = ^^{H+-El+) — 



Wg 



+ ijG'{H--EL-) 



1 



(89) 



2 dR+i,dR- 



(90) 



+max(Ei+,Ei-) - 
+ ^[El^ - max{EL 

+ ^[El~ - max{EL 



(91) 



+ 



+ 



i?L-)]e(^^''-^")-^"- 
i?L-)]e(^^"-^^)-^^- 



(92) 



(93) 
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where the operators H"^ are both identical to H , except that and H act on the space 
of positive and negative configurations, respectively. 

^ H{n^) = -\Y.^ + ^(R^)- (94) 

The coefficients c^,y are real coefficients and will be defined below. To justify that this 
operator is the diffusion operator corresponding to FMC with no cancellation process, we 
need to check that the short-time dynamics described by 



iD-ET)Iit (95) 



dt 

is indeed realized via the two first steps of the FMC algorithm (Langevin and branching 
steps). In the expression of D, the two operators appearing in Eq. (l89l) define a Fokker 
Planck operator in analogy to Eg. (1401) . This operator is the diffusion operator associated 
with the Langevin process, Eq. fl64l) . The term in Eq. flQOj) is a coupling term between the 
moves of positive and negative walkers taking into account the correlation of the gaussian 
random variables rj'^ and rj~. The quantities c^,^(R"^, R") introduced in Eq. flQOj) are nothing 
but the covariance of these variables 

c^,(R+,R-) = (r/+^;). (96) 

The three last contributions describe the branching processes at work in FMC. One recog- 
nizes in Eq.f lQll) the branching of a pair, Eq.( l69i) . The two following contributions, Eqs.( l92|l 
and (l93ll . correspond to the creation of pairs (R^,PR^) and (PR^,R^), with the respec- 
tive weights given by ([70]) and Note that in Eqs. (l92]l93D the operator e^^^''^^^-^^ is 
written in a symbolic form representing a translation of the vector R+ to PR~, the action 
of this operator on the pair (R^, R^) being indeed to create the pair (PR~, R^). 

Now, let us prove that the lowest eigenvalue of D is Eq (bosonic ground-state). For 
that purpose, it is convenient to introduce the operator TZ which transforms a distribution 
of pairs of walkers into a distribution of walkers and then to define the following reduced 
density 



imAK] = / dR' 



UtiR^R') ^ nt(R',PR) 



V^S(R) ■ ^c(R) ■ ^^^^ 
The density TZHt represents the sum of the distributions sampled by each type of walkers 
when the contribution of the other type of walkers is integrated out. Using the explicit 
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expression of D it is a simple matter of algebra to verify that 



TZDUt = HTZUt. (98) 
Using Eqs. (!98|) and (p5|) . one can also write 

~{H-ET)nUt (99) 



dnUt 



dt 

which means that the reduced density evolves under the dynamics of H. In other words, the 
set of positive and negative walkers sample the same distribution as a standard Diffusion 
Monte Carlo algorithm. Now, suppose that Xs is the lowest eigenvalue of D and the 
corresponding eigenstate (the stationary density of the process described by Eq.f P^ ) 

DUs = XsUs. (100) 

Applying 71 on both sides of this identity and using the relation (1^ one gets 

HTZUs = XsTlUs. (101) 

In other words, the reduced density T^fl^ is a positive eigenstate of H with eigenvalue Xs- 
The bosonic state being non-degenerate, we can conclude that Xs = Eq . This ends our 
proof. 

Let us now consider the genuine FMC diffusion operator including the cancellation pro- 
cess, DpMc- To simplify the notations let us suppose that we are in the symmetric case for 
which ipQ = TpQ = ipG (the common guiding function is symmetric under permutation of 
particles). This particular case is much simple because when walkers meet, there is a full 
cancellation and no residual branching. From an operatorial point of view the cancellation 
step consists in introducing a projection operator, Pc, at each step of the dynamics 

Pc = [l- J dR\ R-^ = R,R~ = R){R+ = R,R- = R |]. (102) 

where | R'^,R~) denotes the usual tensorial product. The full FMC diffusion operator can 
thus be written as 

Dfmc = PcD. (103) 

It is important to realize that the -Dfmc operator defined via Pc and D (Eqs. (11021) and 
(!93l) ) represents indeed an equivalent operatorial description of the stochastic rules of FMC 
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described in section IV. B (Langevin, branching and cancellation steps). Note also that using 
the expression (I103P of -DpMC, we have a simple alternative way of recovering the proof just 
presented above that FMC is a bias-free approach. Since this is an important point of 
this work, let us present this alternative proof. The action of the projection operator Pf. 
is to remove from the sample components of the form | = R, R~ = R) for which the 
antisymmetric component of the reduced density is zero 

An\R+ = R, R- = R) = ^— i— (I R)+ I PR)) = 0. (104) 

In other words one has the following algebraic identity 

ATZPc = An. (105) 

In the general case where ipQ ^ ipQ, the cancellation procedure still corresponds to define 
a new operator written as in Eg. (11031) with P^ satisfying the same identity as in (11051) . 
Applying ATZ to the L.H.S. and R.H.S. of equation fl8^ . one has 

dATZUt 



dt 



-ATZOpMc^t = -A{H - ET)nUt. (106) 



This equation indicates that the evolutions of the antisymmetric component of the reduced 
density under the dynamics of -Dfmc and H are identical. This confirms that the energy 
estimator, Eq.( |571) . or any observable estimator not coupling directly positive and negative 
walkers, is not biased. In the case of the energy, the estimator can be written as a function 
of the reduced density as follows 

Let us now turn back to our discussion of the stability of FMC. For that, we need to 
compare the lowest eigenvalue of -Dfmc and, , the lowest eigenvalue of H. Now, it is 
clear from the definition of D^., Eq. (ll03p . that the following relation holds 

Ei > E^. (108) 

Indeed, the action of Pc present in the definition of D^, Eq. (ll03p . consists in removing 
positive coefficients within the extradiagonal part of Z^. As well-known, a consequence of 
such a matrix (or operator) manipulation is to increase the energy of the lowest eigenvalue of 
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the matrix. Expressed in a more physical way, the cancellation process reduces the growth 
of the population of pairs: e~^^o-^T)t ^ ^-{eB-EtY _ 

To summarize, we have shown that FMC reduces the instability of fermion simulations. 



FMC operator, Z^fmc- Because of the inequality fllOSp . this ratio decreases slower in FMC 
that in any standard transient DMC or nodal release methods. We have shown that the 
cancellation process is at the origin of this improvement; however, as we shall see in the 
next section, the cancellation process is efficient (i.e., we have a small difference — Eq) 
only if the correlation between walkers described by the coupling terms c^jy is introduced. 
This feature is important, particularly in high- dimensional spaces where the probability 
of meeting and cancelling becomes extremely small for independent walkers. As a result, 
the correlation of positive and negative walkers is a fundamental feature of FMC. The 
quantitative effect of the correlation on the stabilization of the algorithm is not easy to study 
theoretically and to optimize in the general case. In the next section we will give a numerical 
illustration, for a simple system, of the interplay between cancellation and correlation (via 
the c^i, parameters), and also of the role of the choice of the guiding functions, ipc V^g- 

VI. NUMERICAL STUDY 

A. The model: 2D-harmonic oscillator on a finite grid 

In this section we study the FMC method on a very simple model on a lattice. For this 
model it is possible to calculate Eq (fermionic ground-state energy), Eq (bosonic ground- 
state energy), and Eq the lowest eigenvalue of the FMC operator by a standard deterministic 
method (exact diagonalization). The results obtained for this simple model will provide us 
with a well-grounded framework to interpret the Fermion Monte Carlo simulations. The 
second motivation is that, using such a simple model, it is possible to study the limit of 
large number of walkers, large with respect to the dimension of the Hilbert space considered. 
This possibility turns out to be essential to better understand the FMC algorithm. 

Our model is based on the discretization of a system describing two-coupled harmonic 
oscillators 



The signal-over-noise ratio decreases as e^^o 



where Eq is the lowest eigenvalue of the 




(109) 
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with 



V{x,y)^]^x' + ]^\y'' + xy (110) 



In the foUowing we shah take A = 2. Now, we define the discretization of this model on a 
NyiN regular grid {N odd). A grid point i e (1, . . . , N"^) has the following coordinates 

N , , N 



R 

where 



+ 1)5,, (-- + /- 1)5,) ke[l...Nl le[l...N] (111) 



5x = 5y = Xmax/N (112) 

On this lattice the Hamiltonian has a corresponding discrete representation given by a finite 
matrix. The diagonal part of the matrix reads 

Hu^^ + ^2+y{Ri) ^e(l,■■■,A^') (113) 

Ox Oy 

and the off-diagonal part reads 

^ij — FT when Ri and Rj are nearest-neighbors on the lattice 

Hij = otherwise (114) 

This Hamiltonian is symmetric with respect to the inversion P of center O — (0, 0). 

P{x,y) = {-x,-y) (115) 

As a consequence, the eigenfunctions are either symmetric or antisymmetric under P. We are 
interested in the energy Eq of the lowest antisymmetric eigenstate, 4>q. Even for this simple 
system, we are confronted with a sign instability and a genuine "sign problem". Indeed, 
the sign of (pQ for each grid point cannot be entirely determined by symmetry. Symmetry 
implies only that 0^ vanishes at the inversion center and that the two-dimensional pattern 
of positive and negative values for ([Jq is symmetric by inversion. The precise delimitation 
between positive and negative zones of the wavefunction (analogous to "nodal surfaces" for 
continuous systems) is not known. 

Let us now introduce the trial functions ipr and ips- V'r has to be an (antisymmetric) 
approximation of 0^ and ips some symmetric and positive approximation of the lowest 
eigenstate, 0f . We have chosen them as discretizations of the exact solutions of the initial 
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continuous model. To find these solutions we perform a diagonalization of the quadratic 
form of the potential 



11 11 

y) = ^ + xy^ -kix^ + 2^2f- (116) 



It is trivial to verify that 



^ _ cos^ 9 - X sin^ 6 
^ " cos2e 

A cos^ 6 - sin^ 



tan 29 = 



cos29 
2 



A- 1 



with 



and 



X = xcos^ — ysin6' (117) 



y = xsin^ + ycos^. (118) 

If ki < k2 we choose as trial wavefunction 

Vk2~2\ 



i^T = X exp {-^x^ - ^f), (119) 



while, in the other case, we take 



kl ^2 Vk2~2\ 



= ^exp {-^x' - ^y') (120) 



The lowest (symmetric) eigenstate is chosen to be 



ki ^2 Vk2~2\ 



= exp(-^£^ - ^y'). (121) 

Note that, in the limit of a very large system the trial functions, ■^r and %l)s reduce to two 
exact eigenstates of H; however, this is not the case for finite systems. 

B. FMC on the lattice 

Before presenting our results, let us say a few words about the implementation of the 
FMC on the lattice. The same ingredients as in the continuum case hold, except that in the 
lattice case the Langevin process is realized through a discrete transition probability matrix. 
The probability for a (positive or negative) walker i to go to j, after a time step r is 
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where r is small enough to have a positive density, namely 



1 . ^ 

r < . (123) 



The local energies are defined as in the continuum, Eq.(H3 



• ± 



with the same expression for the guiding functions ipa^ Eq. (H7|) 



E^K) = -^(R) (124) 

Wg 



V'g(R) = V# + ± ciPt. (125) 

Let («i, ^2) represents a given pair of positive and negative walkers (Rj, Rjj)- ^ standard 
diffusion Monte Carlo (no correlation and no cancellation of positive and negative walkers), 
the density of pairs of positive and negative walkers Tl^^l^ evolves as follows in one time-step 



= E n£,P+(ji - ^l)W;^P-iJ2 - ^2)Wr^ (126) 



hh 

where is the Feynman-Kac weight, Eq.( l45l) . To build the FMC algorithm, one first 
correlate the two stochastic processes and P~, Eq.f ll22p . The way it is performed here 
is the counterpart of the correlation term introduced by Liu et ai. 20] in the continuum 
case, Eq. (!63|) . With such a choice, the positive and negative walkers of a pair tend to get 
closer or to move away in a concerted way. For the lattice case it is done as follows. The 
positive walker ji is connected by to a finite number of states 2\ (here, maximum five) 
with probability P^{ji jf). The negative walker j2 is connected to a finite number of 
states J2 with the probability P^{j2 — ^ ji)- The states are ordered taking as criterium 
the distance to the negative walker j2, |Rjc — Rjjl- We do the same for the states ^2 
ordered by their distance with respect to the positive walker. An unique random number 
uniformly distributed between and 1 is then drawn and the repartition functions of the 
two probability measures, p{ji) = P^{ji ji) and ^(^2) = P^ij2 32) ^'^^ then sampled 
using this common random number. The new pair {jlij^ is drawn accordingly. Such a 
procedure defines a correlated transition probability in the space of pairs 

Pc{jij2 ^ ^1^2) 7^ P{ji ^ ^2)P{Jl ^ ^2) (127) 

whose role is to enhance the probability of having positive and negative walkers meeting at 
the same site. Note that, by construction, the correlation introduced via Pc does not change 
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the individual (reduced) densities associated with each type of walker (positive/negative). 
Now, let us write down explicitly the FMC rules in our lattice case, that is, the one time-step 
{k ^ k + 1) evolution of the density of pairs ^jlj^'- 

(i) Correlation and branching. The branching of an individual pair {ji,j2), Eq. fpTl) . 
corresponds to the following evolution of the density 



n 



(fc+i) 

1112 



r(fc) 



(128) 



JU2 



The creation of pairs, Eqs. (l92|93p . can be written as follows 

= Kkr) + E ni'iPc(jij2 - ^l^2)e[wt„ w-]\wt, - w-\/2 



n 



(fc+i) 

-P(«2)«2 



JU2 



(129) 
(130) 



where 0{x,y) = 1 if x > y, 9{x,y) = 0, otherwise. 

(ii) Cancellation. The cancellation process is done when a positive walker and a negative 
walker meet ii = 12 = i 



^(fc+i) 



l-min(^(z),%(z)) 



G 



:i3ii 



If V'g(^) > V'g(^) have 



, [1 - V^G(0/^G(0] -rTW 



(132) 



If ipcii) < 'ipci'^) have 



n 



(fc+i) _ Tr{fc) , [1 - ^g(0/^g(0] 



iP(i) 



(133) 



The operations fll28|129lll30lll31|132lll33l) describe the one time-step dynamics of the sim- 
ulation. At iteration k, the distribution of pairs 1,^2) is obtained and the transient 
estimator of the energy ( !58l) can be computed from 



(fc) imprih) _ gyrfe) ] 



^«l,«2 -'^«l,J2L^(j^) Vg(«2) 



(134) 



This estimator converges to Eq when k —* 00. 

Now, it is important to realize that the FMC rules just presented have, in principle, 
no stochastic character at all. For a finite system the FMC rules can be viewed as simple 
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deterministic matrix manipulations between finite vectors of the Hilbert space (liere, tlie 
multiplications to be performed have been explicitly written). At the beginning of the sim- 
ulation (iteration k = 0) some arbitrary starting vector Tlf^\^ is chosen and, then, iterations 
are performed up to convergence. This important remark will allow us to organize our dis- 
cussion of FMC results into two parts. In a first part (Section VI. C), we perform explicitly 
the matrix multiplications involved and any stochastic aspect is removed from the results. 
Stated differently, this procedure can be viewed as performing a standard FMC simulation 
with an infinite number of walkers (the distribution at each point of the configuration space 
is exactly obtained, no statistical fluctuations are present). In a second part (Section VI. D), 
we implement the usual stochastic interpretation of the FMC rules using a finite number of 
walkers. This second part will allow to discuss the important consequences of the finiteness 
of the number of walkers and, in particular, the role played by the use of population control 
techniques. 

C. FMC using an infinite number of walkers: the deterministic approach 

1. No systematic error: FMC is an exact method 

In this section we verify on our simple example that the FMC rules do not introduce any 
systematic error (bias). The energy expression (11341) has been computed by iterating the 
applications of the elementary operators defined by fll28|129|130|131|132lll33l) . In practice, 
this corresponds to iterate a matrix Gfmc(t)- The distribution 11'^'^+^) at iteration A; + 1 is 
obtained from 11^^^ as follows 

n('=+i)=GFMc(r)n« (135) 

The operator Gfmc(''") has been applied a large number of times on some initial density 
^fij2 component vector, being the linear size of our lattice) and the energy f ll34p 

has been computed at each iteration k. We have checked that the energy converges to the 
exact value, Eq, corresponding the lowest antisymmetric state, with all decimal places. We 
have verified that this is true for several cases corresponding to ranging from 4 to 17. For 
this specific problem these results confirm numerically that FMC is an exact method. 
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2. Meeting time between positive and negative walkers 



Here, we want to illustrate quantitatively the fact that the correlation introduced via the 
probability transition helps greatly to lower the meeting time between positive and negative 
walkers. The influence of the choice of the guiding functions (here, parameter c in Eg. (11251) ) 
on the meeting time is also examined. The meeting time is deflned and evaluated as follows. 
We start with a conflguration consisting of a positive walker located at a corner of the lattice 
and a negative walker located at the opposite corner. The positive and negative walkers are 
moving stochastically with the transition probability deflned in (I122p . We test the two cases 
cases corresponding to uncorrelated and correlated moves. The average time (T) (number of 
Monte Carlo steps times r) before the walkers meet is computed. Our results are presented 
in Table [T] and are given for different linear sizes of the grid. In this flrst case the guiding 
functions are chosen with a large antisymmetric component, c = 4. The results indicate 
clearly that more than one order of magnitude is gained by correlating the moves of the two 
stochastic processes. In Table [III the same calculations are done, except that a symmetric 
guiding function c = {ipQ = i^c — ''Ps) is employed. The average meeting time is found to 
be much lower than in the non-symmetric 4, by nearly two orders of magnitude. 

This is true whether or not the stochastic processes are correlated. This behaviour of the 
meeting time as a function of c is not surprising. When c is large the two functions ipQ and 
ipQ are localized in the nodal pockets of ipT- In the large-c limit ipQ is zero whenever ip^ is 
negative and %pQ is zero whenever %pT is positive. In this limit the overlap between the two 
distributions ipQ and ipQ is zero and we have a similar result for the probability that walkers 
meet. /^From these preliminary results the introduction of non-symmetric wavefunctions 
seems to deteriorate the stability, this property will be conflrmed in the next section. 

3. Stability in time of FMC 

We know from section V that the stability in time is directly related to the magnitude of 
the reduced Bose-Fermi energy gap 

Kb^f = - (136) 

where Eq is the lowest eigenvalue of the FMC diffusion operator. The greater this gap is, the 
faster the signal-over-noise ratio of the Monte Carlo simulation deteriorates, the full stability 
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being obtained only when this gap vanishes. The ultimate goal of an efficient FMC algorithm 
is to reduce the Bose- Fermi gap from its bare value Ab-f = Eq — Eq to a value very close 
to zero (ideally, zero). The energies Eq and Eq can be calculated by exact diagonalization 
of the Fermion Monte Carlo operator, Gfmc- In practice, we have chosen here to extract 
this information from large-time behaviour of the denominator of the energy, Eq. (11341) . In 
this regime the denominator behaves as in Eq. (ITHjl where the reference energy is adjusted 
to keep the number of pairs constant, Et = Eq. 

{V{t = kr)) oc e-(^o^-^<?)* (137) 

The gap Eq — Eq has been extracted from the large-k values of T'(t = A;r), a quantity 
calculated deterministically by iterating the matrix Gfmc- The results for different values 
of c are reported in Table IIIII For both the correlated and uncorrelated processes it is found 
that the gap increases with c. The minimal gap is obtained for c = 0, that is when the 
guiding functions are symmetric ipQ = ipQ = ips- This result is easily explained from the 
fact that there are two factors which favour the cancellation of walkers. First, as we have 
already seen, the average meeting time is minimal when c = since, in this case, the overlap 
between the functions ipQ and tpQ is maximal. Furthermore, a full cancellation between the 
walkers is precisely obtained when c = 0. In conclusion, the greatest stability is obtained 
for a symmetric guiding function. 

In Table [IV] the gaps obtained at c = for different linear sizes are reported. This 
table shows that, in the limit of large grids, that is for a system close to the continuous 
model, the FMC algorithm reduces the gap by a factor ~ 20. Such a result corresponds 
to a huge gain in the stability since projection times about twenty times larger than in a 
standard nodal release method can be used. 

D. FMC using a finite number of walkers: the stochastic approach 

In the previous section the Fermion Monte Carlo method has been discussed and imple- 
mented by manipulating the exact Fermion Monte Carlo diffusion operator without making 
reference to any stochastic aspect (as already mentioned it is formally equivalent to use an 
infinite number of walkers). Of course, for non-trivial systems it is not possible to propagate 
exactly the dynamics of the FMC operator. Accordingly, a finite population of walkers is in- 
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troduced and specific stochastic rules allowing to simulate in average the action of the FMC 
operator are defined. Now, the important point is that in practice -like in any DMC method- 

the population 



22 



2^ This 



one does not sample exactly the dynamics of the FMC operator because o: 
control step needed to keep the finite number of walkers roughly constant, 
step introduces a small modification of the sampled diffusion operator which is at the origin 
of a systematic error known as the population control error. For a bosonic system, the error 
on the ground-state energy behaves as (M is the average size of the population) and an 
extrapolation in can be done to obtain the exact energy. For a fermionic system, as we 
shall see below, this behaviour is qualitatively different and, furthermore, depends on the 
guiding function used. To have a precise estimate of the mathematical behaviour of the 
population control error is fundamental since, in practice, it is essential to be able to reach 
the exact fermi result using a reasonable number of walkers. As we shall see later, this will 
not be in general possible with FMC. 

In this section the Fermion Monte Carlo simulations are performed using 
Eqs.f ll28|129lll30lll31|132p which allow to propagate stochastically a population of M walk- 



ers. The population is kept constant during t 



figuration Monte Carlo (SRMC) method. 
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le simulation by using the stochastic recon- 



2J] In short, the SRMC method is a DMC 



method in which a reconfiguration step replaces the branching step. A configuration step 
consists in drawing M new walkers among the M previous ones according to their respective 
Feynman-Kac weight (for the details, see the references given above). 
In Figure [H the time-averaged energy defined as 

^ Ifl 

is ploted as a function of K. In this formula M{k) and 'D{k) represent the numerator and 
denominator at iteration k of the estimator (I75]l evaluated as an average over the population 
of pairs. In Figure [2] we plot the time-averaged denominator given by 

VK = ^jZ'^{k). (139) 
^ k=i 

The time dependence of this quantity is interesting since it can be used as a measure of 
the stability of the algorithm. |l8| As we have shown above, the algorithm is stable only 
when the reduced Bose- Fermi energy gap, ^b-f = Eq — Eq is equal to zero. Equivalently, 
the denominator (I139p must converge to a constant different from zero. In our simulations 
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the number of walkers was chosen to be M = 100, a value which is much larger than the 
total number of states of the system (here, nine states). Of course, such a study is possible 
only for very simple systems. As seen on the Figures [T][2l and [3] the results obtained in 
the case c = (symmetric guiding function) and c 7^ are qualitatively very different. 
Figure [T] shows that, within statistical error bars, there is no systematic error on the energy 
when a symmetric guiding function is used, c = 0. However, the price to pay is that the 
statistical fluctuations are very large. This point can be easily understood by looking at the 
behaviour of the denominator, Figj2l Indeed, this denominator vanishes at large times, thus 
indicating that the simulation is not stable. In sharp contrast, for c = 4 (non-symmetric 
guiding functions), the statistical fluctuations are much more smaller (by a factor of about 
40) but a systematic error appears for the energy. Furthermore, the denominator ploted 
in Figl2] is seen to converge to a flnite value. The stability observed in the case c = 4 
seems to conflrm the results of Kalos et ai.jlS] for non-symmetric guiding functions {c 0). 
However, the situation deserves a closer look. Indeed, the existence of this flnite asymptotic 
value seems to be in contradiction with our theoretical analysis: the denominator should 
converge exponentially fast to zero, and the algorithm should not be stable {Ab-f > 0). In 
fact, as we shall show now, the asymptotic value obtained for c = 4 and the corresponding 
stability result from a population control error. To illustrate this point, we have computed 
the average denominator as a function of the population size M. Results are reported in 
FiglHl On this plot we compare the population dependence of the denominator (1139^ for 
c = 4 and for a much smaller value of c = 0.5. The values of M range from M = 100 to 
M = 6400. A flrst remark is that the population control error can be quite large and is 
much larger for c = 4 than for c = 0.5. In the Appendix it is shown that the theoretical 
asymptotic behaviour of the error as a function of M is expected to be in 1/M. In the 
c = 0.5 case, the denominator is clearly seen to extrapolate to zero like In the c = 4 case, 
we can just say that the data are compatible with such a behaviour but even for the largest 
M reported in Figj3] (M = 6400) this asymptotic regime is not yet reached. Much larger 
populations would be necessary. This result illustrates the great difficulty in reaching the 
asymptotic regime, even for such a simple system having only nine states. Stated differently, 
the stability observed when using non-symmetric guiding functions disappears for a large 
number of walkers, thus confirming that the stability obtained at finite M is a control 
population artefact. Note that a large population control error on the denominator is not 
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surprising. Indeed, when c 7^ 0, the local energies of the guiding functions, Eq. (H3!) . have 
strong fluctuations because iI}q is far from any eigenstate of H {ipQ contains a symmetric 
and an antisymmetric components). In the case of a symmetric guiding function (c = 0), 
the distribution of walkers is also symmetric at large times and, thus, the average of this 
distribution on the antisymmetric wavefunction ipf must necessarily be zero. Consequently, 
in the c = case there is no control population error on the denominator, FiglJl 

An example of the behaviour of the energy as a function of the population size M is 
presented in figure HI Some theoretical estimates of the energy bias dependence on M are 
derived in the Appendix. Let us summarize the results obtained. When c = (use of a 
symmetric guiding function), the systematic error behaves as in a standard DMC calculation 
for a bosonic system 

5E oc ^7- (140) 

However, the statistical fluctuations are exponentially large since the calculation is no longer 
stable. Now, when c > the systematic error has a radically different behaviour. For a 
population size M the control population error grows exponentially as a function of the 
projection time t 

(X ^e*^--^ (141) 

where I^b-f is the reduced Bose- Fermi energy gap. This dependence of the control pop- 
ulation error as a function of the projection time is of course pathological and is a direct 
consequence of the use of non-symmetric guiding functions. Now, because of the form ( I14ip 
it is clear that a population size exponentially larger than the projection time is necessary to 
remove the systematic population control error. In practical calculations, for a given popu- 
lation of walkers M, one has to choose a finite projection time, t. This time has to be small 
enough to have a small finite population control error but, at the same time, large enough 
to extract the exact fermionic groundstate from the initial distribution of walkers. The best 
compromise is easily calculated and leads to the following expression of the systematic error 
as a function of the number of walkers (for more details, see the Appendix) 

oc — (142) 

where 

7 = ^ — — < 1 (143) 
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and Ai? is the usual Fermi gap (energy difference between the two lowest fermionic states). 

In figure H] some numerical results for the c = 4 and c = 0.5 cases are presented. The 
number of walkers considered are M = 100, 200, 400, 800, and M = 1600. No data are shown 
for the symmetric case, c = 0, because of the very large error bars, see FigUl The calculations 
have been done for the smallest system, = 3 (recall that the finite configuration space 
consists of only nine states) and for very large numbers of Monte Carlo steps (more than 
10®). As it should be, the systematic errors are found to be larger for the c = 4 case than for 
the c = 0.5 case (note that the data corresponding to M = 800 and M = 1600 must not be 
considered because of their large statistical noise). The concavity of both curves confirms 
our theoretical result, 7 < 1, Eqs. fll42fl43p . However, it is clear that to get a quantitative 
estimate of this exponent is hopeless because of the rapid increase of error bars as a function 
of M. The only qualitative conclusion which can be drawn by looking at the curves is that 
7c=4 < 7c=o.5) in agreement with our formula, Eq. fll43p . Finally, let us insist on the fact that, 
despite these very intensive calculations for a nine-state configuration space, no controlled 
extrapolation to the exact energy is possible. 

To summarize, when ipc has an antisymmetric component, the error is expected to de- 
crease -for M large enough- very slowly as a function of the population size [algebraically 
with a (very) small exponent], while in the symmetric case the bias has a much more in- 
teresting -^-behaviour. However, in this latter case the price to pay is the presence of an 
exponential growth of the statistical error. In both cases, and this is the fundamental point, 
the number of walkers needed to get a given accuracy grows pathologically. In addition, 
as illustrated by our data for the very simple model problem treated here, the asymtotic 
regimes corresponding to the l/M'^'-behaviour appear to be very difficult to reach in practice 
(very large values of M are needed). 



VII. CONCLUSION AND PERSPECTIVES 



The FMC method differs from the DMC method by correlating the diffusion of the walkers 
and introducing a cancellation procedure between positive and negative walkers whenever 
they meet. In this work we have shown that the Fermion Monte Carlo approach is exact 
but in general not stable. FMC can be viewed as belonging to the class of transient DMC 
methods, the most famous one being probably the nodal release approach. jQ, [3^ However, 
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in contrast with the standard transient methods, FMC allows to reduce in a systematic way 
the fermi instability. The importance of this instability is directly related to the magnitude 
of some "effective" Bose- Fermi energy gap, A^^^ = — , where E^ is the exact 
Fermi energy and Eq some effective Bose energy. We have seen that this gap is intimately 
connected to the the cancellation rate, that is to say, to the speed at which positive and 
negative walkers cancel. We have shown that E^ < Eq < Eq , where Eq is the standard 
bosonic ground-state energy. As an important consequence, the closest Eq is from the 
exact fermionic energy, the smoother the sign problem is. For the toy model considered, 
the reduction obtained for the instability is very large (orders of magnitude). For large 
dimensional systems, there are strong indications in favor also of an important reduction. 
A first argument is purely theoretical. In FMC the walkers within a pair (R^,R~) are 
correlated in a such a way that — R~ makes a random move only in one dimension. As a 
result there is a high probability that the walkers meet in a finite time even if they move in a 
high-dimensional space. The second argument is numerical. As shown by previous authors, 
the impact of correlating walkers on the average meeting time is important even for much 



larger systems. 18l. l21| 



We have shown that the recent introduction of nonsymmetric guiding functions in FMC 
introduces a large systematic error which goes to zero very slowly as a function of the pop- 
ulation size [~ l/Af^, 7 = Ap/{Ab-f + ^f) and Ap = E[ — Eq is the usual fermionic 
gap]. For an infinite number of walkers, this systematic error is removed and the algorithm 
recovers the Fermi instability. Morever, we have shown that using such guiding functions 
does not in general improve the stability. For a large enough number of walkers, the simu- 
lation can be less stable than the simulation using a symmetric guiding function. Finally, 
it is important to emphasize that the conclusion of this work is that the FMC algorithm 
is not a solution to the sign problem. However, it is a promising way toward "improved" 
transient methods. As a transient method, FMC is expected to converge much better than 
a standard nodal release method. We are presently working in this direction. 
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APPENDIX: POPULATION CONTROL ERROR IN FMC 

FMC, like any Monte Carlo method using a branching process, suffers from a so-called 
population control bias. This systematic error appears because the branching rules (cre- 
ation/annihilation of walkers) are implemented using a population consisting of a large but 
finite number of walkers. Nothing preventing the population size from implosing or ex- 
ploding, a population control step is required to keep the average number of walkers finite. 
A standard strategy to cope with this difficulty consists in introducing a time- dependent 
reference energy whose effect is to slightly modify the elementary weights of each walker 
by a common multiplicative factor (close to one) so that the total weight of the population 
remains nearly constant during the simulation. This step, which introduces some correlation 
between walkers and, therefore, slightly modifies the stationary density, must be performed 
very smoothly to keep the population control error as small as possible. In practice, for stan- 
dard DMC calculations done with accurate trial wavefunctions and population sizes large 
enough, the error is found to be very small, in general much smaller than the statistical 
fluctuations. As a consequence, the presence of a population control bias is usually not con- 
sidered as critical. Here, the situation is rather different. In FMC the use of bosonic-type 
guiding functions introduces very large fluctuations of the local energy and the cancellation 
rules a very small signal-over-noise ratio for fermionic properties. In this case, it is not clear 
whether the bias can be kept small with a reasonable number of walkers. 

In this section we present an estimate of the population control bias in FMC. As we 
shall see our estimate shows that the sign problem is actually not solved but attenuated in 
FMC (an exponentially large number of walkers is needed to maintain a constant bias as the 
number of electrons is increased). The derivation presented in this section is very general: 
it is valid for any exact fermion QMC method based on the use of a nodeless bosonic-type 
reference process and some projection to extract the Fermi ground-state. Accordingly, we 
have chosen not to use the specific framework and notations of FMC but, instead, notations 
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of a general DMC algorithm (transient method). Of course, we do not need FMC to introduce 
nonsymmetric guiding functions. The adaptation of what follows to FMC is straightforward. 

In quantum Monte Carlo we evaluate stochastically the following expression for the lowest 
eigenstate of energy Eq 



where /o is some positive initial distribution and tpT an approximation of the eigenstate 
with energy E^. Here, we deal with a fermionic problem so ipT must be antisymmetric. The 
expression (11441) gives the exact energy E^ only when taking the limit t — > cxd. In practice 
for a finite t, there is a systematic error 

AE'p = E'p- < oc e-^^* = e-(^"-^o")* (145) 

where E^ is the energy of the first excited state in the antisymmetric sector and, A^?, the 
fermionic energy gap. For an exact algorithm with one walker (e.g. Pure Diffusion Monte 



Carlo, 
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371]) one computes the R.H.S of (11441) by evaluating the following expression 



H4>T [P,{t)\c- ds{EL\R{s)]-ET)^ 



Ep = ; r- (146) 



t/jC I 

where ipc is the guiding function, strictly positive, with eventually an antisymmetric com- 
ponent. We have also introduced in this expression, the local energy of the guiding function 

E, = f^ (147) 

The integral in (I146p is done over the drifted random walks going from to t. To simplify 
the notations we note (I146P as follows 

where 

h' ^ ^m)] (149) 

W' = e-/o'^^^^[^(^)] (150) 

P' = ^[R{t)] (151) 
Vg 



E'p = (148) 
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For M independant walkers Ri one has 



M "'I 1 



= Y.z: ' :( (152) 



In the analysis presented here based on a population of M the walkers branched according 
to their relative multiplicities, (the M walkers are therefore no longer independant), one 
replace the individual weight Wi by a global weight 2S[ 

W' = j^T.W! (153) 



As a result the energy may be written as 



Ep = ) - ( (154) 



where 



(155) 



M 

P'^^HpI (156) 



Expression fll54p is exact when the weights are included. A control population error 
arises when one does not take into account the weights in expression fll54p . This population 
control error is thus given by 

^< = 1^- (157) 

(158) 

or, after normalizing W* in such a way that (W^) = 1 (for example, by suitably adjusting 
the reference energy Et), 











if) 










cov(/i*,iy*) + 








cov(p*, W^) + 







This is our basic formula expressing the systematic error at time t resulting from the use 
of a finite population. Now, let us evaluate this expression in the large time t and large M 
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regimes. First, we consider the two denominators appearing in the R.H.S. of Eq. fll59p . Let 
us begin with the denominator of the second ratio: 



De = coy{p\ W')+ <p' > (160) 

Because this denominator is nothing but the denominator of the R.H.S. of Eq. fll46p we can 
conclude that D^, does not depend on the population size M and that it vanishes exponen- 
tially fast 

= i^ee-^^-^-^^)*, (161) 

where Et is the reference energy, Et = Eq for a nodal release-type method, and Et = 
Eq > Eq for the FMC method. Let us now look at the other denominator of Eq. fll59p 

Da = {p'). (162) 

This denominator is the usual quantity evaluated during the simulation. It is an approximate 
quantity since it does not include the corrective weights. The asymptotic behaviour of Da 
depends on ipc- We distinguish two cases: 

(i) If ipG is symmetric (c = 0), the stationary density [t large enough) is symmetric. 
Consequently, Da, which is the average of an antisymmetric function, converges to zero 
exponentially fast at large times. For M large enough, in a regime where the dynamics is 
close to the exact dynamics of the Hamiltonian, we know that the convergence is given by 

Da = Ka{M)e-^^-^\ (163) 

where the coefficient Ka depends on M in general. This coefficient will be determined later. 
^From equations fll63p and (11610 . one can evaluate the error on the denominator 

Da-D, = {Ka{M) - Ke)e-^^-^K (164) 

We also know from the definitions of Da and Dc[ (ll62p . fll60p ] that the difference Da — D^ is 
a covariance of two averages 

Da-De = cov(p*, W') (165) 
which, due to the central-limit theorem, behaves as 

Da-D, = ^C(t) (166) 
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where C{t) is some function of t. Identifying fll64p and (11661) . one finally obtains a deter- 
mination of Ka- Finally, we obtain the following behaviour for the systematic error on the 
denominator 

Da-D, = cov(p*, W) cx (167) 

(ii) If ifjQ is not symmetric (c 7^ 0) the stationary density has an antisymmetric component 
and Da converges to a constant different from zero at large times. Of course, this constant 
depends on the number of walkers M. This dependence can be easily found by using the 
central limit theorem as before 

Da-D, = coy{p\ W) = K^. (168) 

Finally, we have just proved that, when the guiding function is not symmetric, the asymp- 
totic behaviour of the denominator is Da oc j^. This important result is in agreement with 
the numerical data shown in figure (j3]). Using exactly the same arguments, the asymp- 
totic behaviour (large M, large t) of the difference of the two numerators in the R.H.S. of 
expression (11591) is found to be the same as Da — -De 

Na-N, = COY{h\ W) ^Da- D,. (169) 

Now, we are ready to write down the expression of the systematic population control 
error, AE|;^ For M large enough, AE^ is well approximated by its first-order contribution 
in the expansion. Here also, we need to distinguish between the nature of the guiding 
function 

(i) If ipG is symmetric one easily obtains 

A< oc ^ (170) 

(ii) If ipG has an antisymmetric component 

A^A^ oc j^e^"-^' (171) 

Let us now evaluate the total systematic error resulting from using a finite time t and a 
finite population size M 

AE^ = AE^^ + AE'p (172) 

where AEp, the systematic error coming from a finite simulation time, is given by Eq. (ll45p 
and AEp is the error just discussed. The strategy consists in determining, for a given 
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systematic error AE^ ~ e, what is the time t and the number of walkers M one should 
consider. The condition for the total systematic error to be of order e is that both terms in 
(11721) are also of order e 

AE'p ~ e (173) 
AE^' ~ e. (174) 

This is true because no error compensation are present, AEp and AEp being generally of 
the same sign (both positive). Our numerical results on the toy model give an illustration 
of this property. /^From both equations (1174^ and (I145P one can deduce the simulation time 
corresponding to such a systematic error 

(175) 

In other words, to obtain an error of order e it is sufficient to stop the simulation at a time 
t of order fll75p . Now, let us come to the number of walkers M needed. If tpc is symmetric, 
we already know from expression (11701) that the systematic error does not depend on the 
projection time and that the number of walkers M and the systematic error e are related as 
follows 

M oc ^. (176) 

If ipG is not symmetric, the equation (I174p can be easily solved. Replacing in Eq. (11741) . 
AEp by its expression (I17ip and using the relation (11751) one finally finds the number of 
walkers required to obtain a systematic error e. 

M oc e~^^~\ (177) 

Let us make some important comments. First, note that in this formula the dependence 
on the guiding function is not included in the exponent, only in the prefactor. Second, 
this formula shows that the FMC algorithm reduces the systematic error by lowering the 
exponent. As already mentioned, the gap is equal to Eq — Eq in a standard release node 
method and equal to Eq — Eq < Eq — Eq in FMC. Third, in the zero-limit gap, the 
behaviour of the systematic error is recovered. This formula shows that the number of 
walkers needed for a given accuracy, e, grows exponentially with respect to the number of 
electrons. Indeed, although the gap is indeed reduced by FMC, there is no reason not to 
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believe that it will still be proportional to the number of electrons. In consequence, the 
"sign problem" fully remains in FMC 

Finally, let us write the systematic error as a function of the finite population M by 
inverting the preceding equation fll77p 

e oc M-^ (178) 

where 

7 = ^ — — (179) 

This latter equation shows very clearly the respective role played by the Fermi gap, A^, and 
the reduced Bose- Fermi gap, A^-f- 
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FIG. 1: A'^ = 3 Energy estimator as a function of the projection time (number of iterations K). 
Comparaison between c = (large fluctuations) and c = 4 (small fluctuations). The exact energy 
is = 1.86822... Number of walkers M = 100. Number of Monte Carlo steps: 4.10'''. 
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FIG. 2: N = 3 Denominator as a function of the projection time (number of iterations K). 
Comparaison between c = 4 (upper curve) and c = (lower curve). Number of walkers M = 100. 
Number of Monte Carlo steps: 4.10^. 
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FIG. 3: A'^ = 3, time-averaged denominator, Eq. (|139p . as a function of 1/M. 
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FIG. 4: = 3 Energy, Eq. (jl34p . as a function of 1/M for the c = and c = 4 cases. Exact energy: 
£^^=1.86822.. .(horizontal sohd hne) 



TABLES 



TABLE L Average meeting times for the correlated and uncorrelated cases in the non-symmetric 
case c = 4, Eq. (jl25p . In this example, Xmax = 3, Eq. ()112p . and r = Q.^Tmax, where Tmax is the 



maximal time-step defined in Eq. (jl23p . 



Linear size 


^ Uncorr. 


^ Lorr. 


iV = 3 


2152(23) 


134(1) 


iV = 5 


2162(20) 


92(1) 


N = 7 


2234(26) 


75(1) 


iV = 9 


2834(27) 


76(1) 


= 11 


3546(38) 


82(1) 


= 13 


4214(41) 


88(1) 


N = lb 




94(1) 


N = n 




102(1) 



TABLE II: Average meeting times for the correlated and uncorrelated cases in the symmetric 
case (c = 0, symmetric guiding function), Eq . (|125p . In this example, Xmax ~ 3, Eq.([Il2|), and 
r = 0.9xmax) where r^ax is the maximal time-step defined in Eq. (|123p . 



Linear size N 


^ Uncorr. 


^ Uorr. 


N = ?, 


3.32(2) 


1.634(9) 


N = 5 


5.64(6) 


2.27(1) 


N = 7 


7.6(1) 


2.65(1.6) 


N = 9 


10.3(1) 


3.32(2.5) 


N = 11 


12.95(8) 


4.18(4) 


N = 13 


15.7(1) 


4.78(4) 
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= 15 


18.5(2) 


5.52(6) 


= 17 


21.6(2) 


6.26(7) 



TABLE III: A'" = 3. Reduced Bose-Fermi gap Ab-f, Eq. (|136p . with or without correlation for 
different values of c. The average meeting times are also indicated in parentheses. The bare 
Bose-Fermi gap, Ab-f, is ~ 0.7695 (here, Xmax = 3. and r = 0.09rmaa;). 



Value of c 


Correlated process 


Uncorrelated process 


c = 
c= 1 
c = 2 
c = 3 
c = 4 


Ab-f = 0.0366 1.634(9)] 
Ab-f = 0.0917 6.37(7)] 
Ab-f = 0.1336 24.6(2)] 
Ab-f = 0.1026 62.9(3)] 
Ab-f = 0.1092 134(1)] 


Ab-f = 0.1629 ; 3.32(2)] 
Ab_^ = 0.2540 ; 16.5(2)] 

Ab_p = 0.2277 ; 130.5(5)] 
AB„f = 0.1981 ; 627(3) ] 

Ab-f = 0.1787 ; 2152(23)] 



TABLE IV: Comparison between the reduced Bose-Fermi gap, Ab-f, and the bare Bose-Fermi 
gap, Ab-f, in the symmetric case (c = 0, symmetric guiding function) as a function of 



Value of N 




Ab- 


-F 




Ab- 


-F 


Gap ratio 


N = 3 


Ab- 


-F = 


0.0366 


Ab- 


-F = 


0.7695 




N = 5 


Ab- 


-F = 


0.0516 


Ab- 


-F = 


1.0195 




N = 7 


Ab- 


-F = 


0.0577 


Ab- 


-F = 


1.1782 


^B-F ^ 0.0490 



49 



